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ABSTRACT 

We  consider  the  problem  of  determining  the  temperature  history  inside  a  gun 
barrel  from  embedded  thermocouple  measurements  at  some  distance  away  from  the 
inside  wall.  This  inverse  problem  leads  to  an  improperly  posed  initial  value 
problem  for  a  nonlinear  system  of  partial  differential  equations,  whenever  the 
thermal  properties  are  temperature  dependent.  We  discuss  a  step-by-step 
marching  algorithm  for  the  numerical  computation  of  such  problems.  The  scheme 
is  stabilized  by  appropriately  filtering  in  the  frequency  domain  at  each  step. 
We  illustrate  this  technique  with  a  numerical  experiment  on  a  nonlinear 
problem  whose  exact  solution  is  known.  The  basic  ideas  are  applicable  to 
other  unstable  evolution  equations^ 


I .  Introduction 

This  report  summarizes  the  results  of  an  important  computational 
experiment  on  a  nonlinear  inverse  heat  conduction  problem  whose 
exact  solution  is  known.  We  consider  the  problem  of  determining 
the  temperature  history  at  the  inside  wall  of  a  gun  barrel,  from 
embedded  thermocouple  measurements  at  various  points  in  the  annular 
metallic  region  between  the  inner  and  outer  radii  of  the  cannon. 

As  the  shell  is  fired,  a  continuous  trace  is  recorded  at  each 
thermocouple,  providing  temperature  as  a  function  of  time  at  the 
corresponding  fixed  spatial  location. 

*  Research  sponsored  by  the  U.S.  Army  Research  Office 
under  MIPR  No.  ARO  63-82 
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The  present  study  centers  around  a  novel  computational  technique 
designed  especially  for  coping  with  the  nonlinear  case  of  tempera¬ 
ture  dependent  thermal  properties.  It  is  a  sequel  to  [1]  where  the 
linear  quarter  plane  problem  with  constant  coefficients,  was  tho¬ 
roughly  analyzed.  As  was  shown  there,  in  that  case,  the  inverse 
problem  can  be  formulated  either  as  a  Volterra  integral  equation 
of  the  first  kind,  or  equivalently,  as  an  initial  value  problem 
for  the  one  dimensional  heat  equation  run  sideways.  Either  formu¬ 
lation  leads  to  an  improperly  posed  problem  in  which  the  solution, 
when  it  exists,  depends  discontinuously  on  the  data. 

n 

The  inverse  problem  can  be  regularized  in  the  L  norm  by  placing  an 

a-priori  bound  M  on  the  norm  of  the  unknown  temperature  history, 

f(t),  at  the  inside  wall  x  =  0;  at  the  same  time,  the  measured  noisy 

temperature  data  g  (t),  at  the  location  x  =  i  >  0,  is  regarded  as 
m 

p 

differing  by  at  most  e.  in  the  L  norm  from  unknown  smooth  exact  data 
g(t),  for  which  a  solution  exists.  It  is  assumed  that  e  and  M  are 
known  and  compatible.  As  shown  in  [1,  equations  (2.20),  (2.21)]  this 
leads  to  explicit  formulae  for  the  temperature  and  gradient  histories 
at  each  fixed  x,  0  <  x  <  i.  Also,  error  estimates  are  obtained  for 
the  regularized  solutions  implying  Holder  continuity  with  respect  to 
the  data,  for  each  fixed  positive  x.  These  estimates  degenerate  at 
the  wall,  [1,  Theorem  1], 

The  regularization  procedure  can  be  interpreted  as  solving  the 
initial  value  problem  for  the  sideways  heat  equation  with  appro¬ 
priately  modified  initial  data.  An  explicit  finite  difference 
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scheme  consistent  with  that  problem  is  shown  to  be  unconditionally 
convergent,  when  used  with  the  filtered  initial  data,  [1,  Theorem  3]. 
This  step-by-step  marching  scheme  in  the  x-variable  is  the  basis  for 
our  approach  to  the  nonlinear  case  of  a  temperature-dependent  dif¬ 
fusion  coefficient.  We  regularize  the  calculation  at  each  step  by 
filtering  in  the  frequency  domain,  using  FFT  algorithms;  we  then 
return  to  the  physical  variables  for  the  calculation  of  the  next 
step.  The  filtering  function  used  at  each  step  is  that  determined 
by  the  related  constant  coefficient  problem.  This  algorithm  is 
outlined  in  [1,  Section  7]. 

In  order  to  test  the  robustness  of  this  procedure,  an  example  was 
manufactured  with  a  known  exact  solution.  This  is  a  fictitious 
mathematical  problem,  artificially  created  so  as  to  have  a  solution 
which  simulates  conditions  presumed  to  exist  in  a  155mm  cannon.  The 
relevant  parameters  were  made  available  to  us  by  Dr.  A.  K.  Celmins, 
U.S.  Army  Ballistic  Research  Laboratory,  Aberdeen  Proving  Grounds, 
Maryland.  The  "exact"  solution  was  constructed  numerically  by 
solving  a  well  posed  direct  problem  as  explained  below. 

2.  The  Direct  Problem 

Consider  the  initial  boundary  value  problem 

3  U  3  3  U 

(2.1)  --  =  -  [a(u)  --],  0  <  x  <  *,  t  >  0, 

at  ax  ax 

(2.2)  u(0,t)  =  f(t),  u( 1 , t )  =  h(t) ,  t  >  0, 
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(2.3)  u(x,0)  =  300°  K 


where  t  is  the  time  measured  in  milliseconds,  x  measured  in 
millimeters  represents  distance  away  from  the  inside  wall,  and 
u(x,t)  is  the  temperature  in  degrees  Kelvin.  The  heat  conduction 
equation  (2.1)  is  a  simplification  of  the  actual  physical  situation, 
in  that  first  order  terms  arising  from  cylindrical  symmetry  have  been 
neglected,  as  well  as  the  variation  of  specific  heat  with  tempera¬ 
ture.  Moreover,  for  gun  steel  at  temperatures  between  300°  K  and 
1000°  K,  the  conduction  coefficient  a(u)  in  (2.1)  is  well  approxi¬ 
mated  by  a  linear  function  of  u, 

(2.4)  a(u)  =  (1.299  -  1.144  x  10-3  (u  .  255)}  x  10-2  mm2/millisec 

We  remark  that  the  methodology  to  be  discussed  can  easily  accommodate 
the  more  exact  differential  equation,  as  well  as  more  complicated 
dependencies  of  a(u)  on  u  .  We  shall  refer  to  the  quantity 

3u 

(2.5)  w(x,t)  =  -a(u)  -- 

3x 

as  the  temperature  gradient,  by  an  abuse  of  terminology.  It  is 
measured  in  mm°  K/mi 1 1 i seconds.  In  all  Figures  shown  below  dealing 
with  plots  of  w(x,t)  as  a  function  of  t  for  some  fixed  x  ,  the 
vertical  axis  bears  the  legend  "temperature  gradient." 

The  functions  f(t)  and  h(t)  in  (2.2)  represent,  respectively,  the 
temperature  histories  at  the  inside  wall  and  at  1mm  away  from  the 
wall.  These  mathematical  functions  are  plotted  in  Figure  1;  they 
are  constructed  so  as  to  approximate  observed  temperature  histories 
in  gun  barrels,  [4]. 
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The  direct  problem  given  by  (2.1),  (2.2),  and  (2.3)  was  solved 
numerically,  using  an  adaptive  partial  differential  equation  software 
package,  M0L1D,  [3].  The  numerical  integration  was  carried  out  to  a 
distance  in  time  equal  to  100  milliseconds.  The  temperature  u(x,t) 
and  gradient  w(x,t)  were  evaluated  at  various  fixed  values  of  x  , 
as  functions  of  time,  and  stored  for  subsequent  comparisons.  Figures 
2  and  3  show  the  histories  of  u  and  w  at  x  =  .25mm.  As  is 
evident  from  Figure  3,  the  numerical  calculation  of  w  is  not  free 
from  noise.  Nevertheless,  we  use  the  term  “exact  solution"  for  any 
history  obtained  by  the  above  numerical  computation  of  the  direct 
problem.  All  histories  are  records  consisting  of  400  equi spaced 
samples  on  the  time  interval  [0,100]  milliseconds. 

3.  The  Inverse  Problem 

The  physical  region  of  interest  here  is  the  x  interval  between 
0  and  .25mm.  The  histories  in  Figures  2  and  3  simulate  what  might 
have  been  recorded  by  a  thermocouple  at  .25mm  away  from  the  inside 
wall  as  the  shell  is  fired.  The  object  is  to  use  such  data  to  re¬ 
construct  the  temperature  and  gradient  histories,  arbitrarily  close 
to  the  inside  wall.  In  actuality,  two  thermocouple  readings  are 
necessary  at  x  =  xQ  and  x  =  Xj,  with  xQ  <  .25  <  x  ;  a  well  posed 
direct  calculation,  as  in  Section  2  above,  then  yields  u  and  w  at 
x  =  .25mm.  As  noted  in  the  references  given  in  [1],  this  type  of 
inverse  problem  occurs  in  a  variety  of  heat  transfer  contexts.  The 
purpose  of  our  computational  experiment  is  as  follows: 

a)  To  demonstrate  the  feasibility  of  the  inverse  calculation 
in  a  realistic  situation  in  which  rapidly  varying  solutions 
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and  nonlinearity  play  a  role.  As  may  be  seen  from  Figure  1, 
the  postulated  temperature  at  the  wall  rises  from  300°  K  to 
almost  1000°  K  in  the  first  10  milliseconds.  In  this 
temperature  range,  the  conduction  coefficient  a(u)  given 
by  (2.4)  undergoes  a  280  percent  change. 

b)  To  demonstrate  the  robustness  of  our  algorithm  with  noisy 
data  and  a  fine  grid. 

c)  The  regularized  marching  procedure  we  shall  use  is  a 
powerful  general  method,  applicable  to  other  ill-posed 
evolutionary  partial  differential  equations,  linear  and 
nonlinear.  As  used  here,  it  is  an  adaptation  to  the  non¬ 
linear  case  of  an  algorithm  which  is  rigorously  justified 
in  the  constant  coefficient  case.  While  the  heuristic 
"local  mode  analysis"  underlying  our  regularization  is 
likely  to  be  valid  in  many  other  cases  of  ill-posed 
initial  value  problems,  there  is  a  need  for  well -documented 
realistic  inverse  calculations. 

Let  z  =  i  -  x  and  let  a  ,  aj  be  positive  constants  such  that 

(3.1)  0  <  aQ  <  a(u)  <  a1  . 

Let  b(u)  =  --  ,  and  let  v(z,t)  =  u(x,t).  Using  (2.5),  we  may  write 
du 

(2.1)  as  an  equivalent  first  order  system 

(3.2)  v  =  -----  ,  w  =  v  ,  0  <  z  <  *,  t  >  0, 

z  a(v)  z  t 

with  the  subscript  notation  used  for  partial  derivatives. 
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to  be  integrated  in  the  direction  of  increasing  z  from  z  =  0  to 
z  =  l;  we  use  the  initial  values  given  in  Figures  2  and  3  and  the 
following  boundary  conditions  at  t  =  0, 

(3-3)  v(z,0)  =  300°  K,  w(z,0)  =  0. 


Let  Az  be  the  increment  in  the  z-variable  and  let  l  =  (N  +  1)az. 
Let  vn(t),  wn(t)  denote,  respectively,  v(nAz,t),  w(nAz,t),  for 
0  <  n  <  N  +  1.  The  following  finite  difference  approximation  is 
second  order  accurate  and  explicit, 

u/jt)  +  “2  vt!t) 
a(*"(t))  2a(vn(t» 


(3.4)  v"+1(t)  =  v"(t)  + 


_  ^Lb(v"(t))Cw"(t)]' 
2a3(vn(t)) 

n  Az2w?(t) 

(3.5)  wn+1(t)  =  wn(t)  +  Azv  (t)  +  - 

1  2a(vn(t)) 


Az2  b(vn(t) )  wn(t)  vj(t) 

2aV(t)) 

An  effective  way  of  implementing  this  scheme  is  to  use  cubic  spline 
interpolation  at  the  400  equally  spaced  mesh  points  on  the  time 
interval  [0,T].  Differentiating  the  spline  function  produces 
0(At  )  accurate  derivatives  v^(t),  w^(t)  at  these  same  mesh  points, 
and  hence  vn+1(t),  wnfl(t)  from  (3.4),  (3.5).  The  next  step  is  to 
stabilize  this  process  by  filtering  each  of  these  functions  in  the 
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frequency  domain.  This  is  accomplished  by  dividing  the  k^h  Fourier 
coefficient  by  the  precomputed  weight  x^,  where 


(3.6)  Xfc  =  (1  +  w2  exp  [  z 


2|k|, 


a0T 


]) 


there  to  =  (-)  is  the  L2  noise  to  signal  ratio.  See  [1]. 
M 


With  z  =  .25mm,  the  x-interval  [0,t]  was  divided  into  450  equally 
spaced  mesh  points,  and  the  above  procedure  was  implemented  with 
to  =  .001.  Figures  4  through  11  summarize  the  comparison  between 
exact  and  computed  solutions  at  the  i nterior  location  x  =  .056mm. 

An  idea  of  the  relative  errors  in  the  calculation  is  easily  gained 
from  Figures  7  and  11.  Although  the  "logarithmic  convexity" 
estimates  in  Theorem  1  of  [1]  degenerate  at  the  wall,  the  computa¬ 
tion  was  pursued  for  450  cycles  and  approximations  to  the  temperature 
and  gradient  histories  at  the  wall  were  obtained.  These  are  shown 
in  Figures  13  and  17.  The  "exact"  temperature  and  gradient  histories 
at  the  wall  are  shown  in  Figures  12  and  16.  Clearly,  slight  in¬ 
accuracies  in  the  well-posed  direct  calculation  of  u(x,t)  near  x  =  0, 
lead  to  a  rather  noisy  determination  of  the  exact  w(x,t)  at  x  =  0;  in 
particular,  the  pronounced  spike  near  t  =  40  milliseconds  in  Figure 
16  is  a  numerical  artifact  which  should  be  disregarded.  Nonetheless, 
we  have  chosen  to  compare  the  computed  gradient  history  in  Figure  17 
with  the  wall  profile  given  in  Figure  16,  As  is  evident  from  Figures 
15  and  19,  the  wall  estimates  obtained  by  solving  the  inverse  problem 
are  quite  reliable.  This  is  especially  true  during  the  first  twenty 
or  so  milliseconds  where  peak  values  are  achieved. 
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4.  Conclusions 


A  regularized  marching  algorithm  has  been  shown  to  be  effective  in 
solving  nonlinear  inverse  heat  transfer  problems  in  gun  barrels.  In 
[23,  a  similar  technique  was  used  successfully  on  linear  backwards 
parabolic  equations  with  highly  variable  coefficients.  More 
recently,  success  has  also  been  achieved  on  other  unstable  examples 
involving  Burgers'  equation  with  the  time  direction  reversed. 

Future  work  should  be  directed  towards  problems  in  two  or  more  space 
dimensions  in  general  domains,  in  the  context  of  heat  transfer  and 
fluid  mechanics. 
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PRESCRIBED  HISTORIES  AT  X=0  AND  X=1 
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FIGURE  1 


TEMPERATURE  HISTORY  AT  .25  NHS 
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FIGURE  2 
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(INITIAL  DATA  FOR  INUERSE  CONFUTATION) 
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FIGURE  A 


TIME  IN  MILLISECONDS 


ESTIMATED  TEMPERATURE  HISTORY  AT  X=.056  MMS 


TIME  IN  MILLISECONDS 
(OBTAINED  BY  SOLOING  INUERSE  PROBLEM) 
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FIGURE  6 


TIME  IN  MILLISECONDS 
COMPARISON  UITH  EXACT  SOLUTION 


71 

O 


_JLU 

-1 

i-«X 


ETf- 


z:3 


I-  UJ  SZ  CL  UJ  (H  <r  h-  3  01 1x1  Q  Ixl  O  01  LU  Ui  CO  V 

A 


FIGURE  7 


ESTIMATED  GRADIENT  HISTORY  AT  X=«056  MIS 
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FIGURE  11 
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FIGURE  12 
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ABSOLUTE  ERROR  IN  ESTIMATED  TEMPERATURE  AT  UALL  (X-0.0) 
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FIGURE  14 


TIME  IN  MILLISECONDS 
COMPARISON  UITH  EXACT  SOLUTION 


GRADIENT  HISTORY  AT  HALL  (X=0«0) 
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FIGURE  19 


TIME  IN  MILLISECONDS 
GRADIENT  IN  MM  DEGREES  K/MILLISECONDS 
COMPARISON  UITH  EXACT  SOLUTION 
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